/**
 * 创建直方图的签名，用于EMD计算
 * 
 */

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;

void help(char** argv) {
    cout << "\nCall is:\n"
    << argv[0] << " modelImage0 testImage1 testImage2 badImage3 \n\n"
    << "for example:"
    << "./ch7_ex7_3_expanded HandIndoorColor.jpg HandOutdoorColor.jpg" 
    << "HandOutdoorSunColor.jpg fruits.jpg\n"
    << "\n";
}

int main (int argc, char** argv) {
    if (argc != 5) {help(argv); return -1;}

    vector<cv::Mat> src(5);
    cv::Mat         temp;
    int             i;

    temp = cv::imread(argv[1], 1);
    if (temp.empty()) {
        cout << "Error on reading image 1 " << argv[1] << endl;
        help(argv);
        return -1;
    }

    // 将第一个图像解析为两个y方向一半大小的
    cv::Size    size        = temp.size();
    int         width       = size.width;
    int         height      = size.height;
    int         halfheight  = height >> 1;

    cout << "Getting Size [[" << temp.rows << "] [" << temp.cols << "]]\n" << endl;
    cout << "Got size (w, h): (" << size.width << "," << size.height << ")" << endl; 

    // 第一种图分为上下两部分，保存在 src[0] 和 src[1]
    src[0] = cv::Mat(cv::Size(width, halfheight), CV_8UC3);
    src[1] = cv::Mat(cv::Size(width, halfheight), CV_8UC3);

    cv::Mat_<cv::Vec3b>::iterator tmpit = temp.begin<cv::Vec3b>();
    
    // top
    cv::Mat_<cv::Vec3b>::iterator s0it = src[0].begin<cv::Vec3b>();
    for (i = 0; i < width*halfheight; ++i, ++tmpit, ++s0it) *s0it = *tmpit;

    // bottom
    cv::Mat_<cv::Vec3b>::iterator s1it = src[1].begin<cv::Vec3b>();
    for (i = 0; i < width*halfheight; ++i, ++tmpit, ++s1it) *s1it = *tmpit;


    // 载入其他3个图
    for (i = 2; i < 5; i++) {
        src[i] = cv::imread(argv[i], 1);
        if (src[i].empty()) {
            cerr << "Error on reading image " << i << ": " << argv[i] << "\n" << endl;
            help(argv);
            return -1;
        }
    }

    // 计算 HSV 图像， 分割为独立的平面
    vector<cv::Mat> hsv(5), hist(5), hist_img(5);
    int         h_bins      = 8;
    int         s_bins      = 8;
    int         hist_size[] = {h_bins, s_bins}, ch[] = {0, 1};
    float       h_ranges[]  = {0, 180};
    float       s_ranges[]  = {0, 255};
    const float* ranges[]   = {h_ranges, s_ranges};
    int         scale       = 10;

    for(i = 0; i < 5; ++i) {
        cv::cvtColor(src[i], hsv[i], cv::COLOR_BGR2HSV);
        cv::calcHist(&hsv[i], 1, ch, cv::noArray(), hist[i], 2, hist_size, ranges, true);
        cv::normalize(hist[i], hist[i], 0, 255, cv::NORM_MINMAX);

        hist_img[i] = cv::Mat::zeros(hist_size[0]*scale, hist_size[1]*scale, CV_8UC3);

        // 画直方图
        for (int h = 0; h < hist_size[0]; h++) {
            for (int s = 0; s < hist_size[1]; s++) {
                float hval = hist[i].at<float>(h,s);
                cv::rectangle(
                    hist_img[i],
                    cv::Rect(h*scale, s*scale, scale, scale),
                    cv::Scalar::all(hval),
                    -1
                );
            }
        }
    }


    // 显示
    cv::imshow("Source0", src[0]);
    cv::imshow("HS Histogram0", hist_img[0]);

    cv::imshow("Source1", src[1]);
    cv::imshow("HS Histogram1", hist_img[1]);

    cv::imshow("Source2", src[2]);
    cv::imshow("HS Histogram2", hist_img[2]);

    cv::imshow("Source3", src[3]);
    cv::imshow("HS Histogram3", hist_img[3]);

    cv::imshow("Source4", src[4]);
    cv::imshow("HS Histogram4", hist_img[4]);

    // 比较直方图
    cout << "Comparison:\n"
        << "Corr            Chi             Intersect       Bhat\n" << endl;

    for (i = 1; i < 5; ++i) {
        cout << "Hist[0] vs Hist[" << i << "]:" << endl;
        for (int j = 0; j < 4; ++j) {
            cout << "method[" << j << "]:" << cv::compareHist(hist[0], hist[i], j) << " ";
        }
        cout << endl;
    }

    // EMD
    vector<cv::Mat> sig(5);
    cout << "\nEMD: " << endl;

    for (i = 0; i < 5; ++i) {
        vector <cv::Vec3f> sigv;
        // 重新正则化直方图，使得总权重为1
        cv::normalize(hist, hist, 1, 0, cv::NORM_L1);
        for(int h = 0; h < h_bins; h++) {
            for(int s = 0; s < s_bins; s++) {
                float bin_val = hist[i].at<float>(h, s);
                if (bin_val != 0) 
                    sigv.push_back(cv::Vec3f(bin_val, (float)h, (float)s));
            }
        }
        // 构造 Nx3 的 32fC1 矩阵，N 是非0直方图bin个数
        sig[i] = cv::Mat(sigv).clone().reshape(1);
        if (i > 0)
            cout << "Hist[0] vs Hist[" << i << "]: "
                << cv::EMD(sig[0], sig[i], cv::DIST_L2) << endl; 
    }


    cv::waitKey(0);
    return 0;
}